%% Calibration
% Runs the calibrations required for each of the models

function [cfg,indepvars, depvars] = Calibration(cfg,exvars,depvars,indepvars,r)
%% Initialize some values before the actual iteration
temp.CalibrationOutput      = [];
modes                       = indepvars.UseModes(r,:,:);

%% Switch between the different models        
switch cfg.Model
    case {'ModelTravelOptimizedTW','ModelTravelLinearTW'}
        %% Run the optimization for TIMER-TRAVEL         
            boundaries.low      = permute([indepvars.MinTimeWeight(r) zeros([1 , sum(indepvars.UseModes(r,1,:))])],[1 3 2]); 
            boundaries.high     = permute([1 10*ones([1 , sum(indepvars.UseModes(r,1,:))])],[1 3 2]);

        if cfg.GeneticSolverT0 % Calibrates the global starting values
            [cfg,depvars,indepvars] = SelectSolver(cfg,exvars,depvars,indepvars,r, cfg.CalibrationYears,{'indepvars.TimeWeight','indepvars.PreferenceFactorStart'},zeros(1,sum(indepvars.UseModes(r,1,:))+1),'ga',boundaries.low,boundaries.high);
        end
        
        %% Refine the Optimization
        [cfg,depvars,indepvars] = ... % Finestune the solution
            SelectSolver(cfg,exvars,depvars,indepvars,r, cfg.CalibrationYears,{'indepvars.TimeWeight','indepvars.PreferenceFactorStart'},...
              [indepvars.TimeWeight(r, cfg.CalibrationYears(1)) permute(indepvars.PreferenceFactorStart(r,1,indepvars.UseModes(r,1,:)),[2,3,1])],cfg.Solver,boundaries.low,boundaries.high);
          
          if cfg.EnableLinearPF==1
              boundaries.PreferenceFactorSlope = -indepvars.PreferenceFactorStart(r,1,:) / length(cfg.TotalTime); 
            [cfg,depvars,indepvars] = ... % Used for the angulation
                SelectSolver(cfg,exvars,depvars,indepvars,r,cfg.CalibrationTime,{'indepvars.PreferenceFactorSlope'},...
                permute(indepvars.PreferenceFactorSlope(r,1,indepvars.UseModes(r,1,:)),[2,3,1]),cfg.Solver,...
                boundaries.PreferenceFactorSlope(indepvars.UseModes(r,1,:))*0, permute([0.1 0.1 0.1 0.1],[1 3 2]));
          end
          
    case {'ModelPoles'}
        %% Version of POLES model based on publication given by Kitous
            boundaries.low         = permute([-0.00001 -0.000000001 -0.0000001],[1 3 2]); 
            boundaries.high        = permute([0.000001 0.000001 0.0000001],[1 3 2]);

         if cfg.GeneticSolverT0
            [cfg,depvars,indepvars]     = SelectSolver(cfg,exvars,depvars,...
                indepvars,r, cfg.CalibrationYears,...
                {'indepvars.TimeTrend'},zeros(1,3),'ga',boundaries.low,boundaries.high);
         end
        
        [cfg,depvars,indepvars]         = SelectSolver(cfg,exvars,depvars,indepvars,r, cfg.CalibrationYears,...
            {'indepvars.TimeTrend'},permute(indepvars.TimeTrend(r,1,modes) ,[2,3,1])...
            ,cfg.Solver,boundaries.low,boundaries.high);
        
    case {'ModelSimple'}
        %% Simplified version of the POLES model, i.e. only elasticities per mode
        

        boundaries.low          = permute([0.2 0.7 0.7 0.7 0.1 0.6 0.6 0.6],[1 3 2]); 
        boundaries.high         = permute([2 2 1.3 2.1 2 1 1 1.3],[1 3 2]);
         if cfg.GeneticSolverT0
            [cfg,depvars,indepvars]     = SelectSolver(cfg,exvars,depvars,indepvars,r, cfg.CalibrationYears,{'indepvars.IncomeElasticity','indepvars.PriceElasticity'},zeros(1,2*sum(indepvars.UseModes(r,1,:))),'ga',boundaries.low);
         end
        
        [cfg,depvars,indepvars]         = SelectSolver(cfg,exvars,depvars,indepvars,r, cfg.CalibrationYears,{'indepvars.IncomeElasticity','indepvars.PriceElasticity'}...
            ,[permute(indepvars.IncomeElasticity(r,1,indepvars.UseModes(r,1,:)) ,[2,3,1]), permute(indepvars.PriceElasticity(r,1,indepvars.UseModes(r,1,:)) ,[2,3,1])],cfg.Solver,boundaries.low,boundaries.high);
        
    case {'ModelGcam'}
        %% Run calibration for GCAM model  
        % First determine sigma over 2005 - this can be done analythically
        % Next calibrate the shareweights
            indepvars.Sigma(r)  = ...
                sum(exvars.TravelDemand(r,cfg.CalibrationYears(end),modes),3)...
                ./((exvars.Income(r,cfg.CalibrationYears(end))^indepvars.IncomeElasticity(r)).* ...
                (sum((exvars.TravelPrice (r,cfg.CalibrationYears(end),modes)...
                + exvars.WageRate (r,cfg.CalibrationYears(end))...
                ./ exvars.Speed(r,cfg.CalibrationYears(end),modes))...
                .*exvars.ModeSplit(r,cfg.CalibrationYears(end),modes),3)^-indepvars.PriceElasticity(r)));
            boundaries.low               = permute([ones([1 ,sum(indepvars.UseModes(r,1,:))])],[1 3 2])*0; 
            boundaries.high              = permute([ones([1 ,sum(indepvars.UseModes(r,1,:))])],[1 3 2])*10000;
         if cfg.GeneticSolverT0 && 0 ==1
            [cfg,depvars,indepvars]     = SelectSolver(cfg,exvars,depvars,indepvars,r, cfg.CalibrationYears,{'indepvars.ShareWeight'},zeros(1,sum(indepvars.UseModes(r,1,:))),'ga',boundaries.low,boundaries.high);
         
        boundaries.low               = permute([ones([1 ,1+sum(indepvars.UseModes(r,1,:))])],[1 3 2])*0; 
        boundaries.high              = permute([ones([1 ,1+sum(indepvars.UseModes(r,1,:))])],[1 3 2])*10000;
        [cfg,depvars,indepvars]         = SelectSolver(cfg,exvars,depvars,indepvars,r, cfg.CalibrationYears,{'indepvars.ShareWeight','indepvars.Sigma'}...
            ,[permute(indepvars.ShareWeight(r,1,indepvars.UseModes(r,1,:)) ,[2,3,1]),permute(indepvars.Sigma(r) ,[2,3,1])],cfg.Solver,boundaries.low,boundaries.high);
         else
             for n =1
             cfg.OFout = false;
             [cfg,depvars,indepvars]         = SelectSolver(cfg,exvars,depvars,indepvars,r, cfg.CalibrationYears,{'indepvars.ShareWeight'}...
            ,permute(indepvars.ShareWeight(r,1,indepvars.UseModes(r,1,:)),[2 3 1]),cfg.Solver,boundaries.low,boundaries.high);
             
            [cfg,depvars,indepvars]         = SelectSolver(cfg,exvars,depvars,indepvars,r, cfg.CalibrationYears,{'indepvars.Sigma'}...
            ,permute(indepvars.Sigma(r,1),[2 3 1]),cfg.Solver,0,1000);
        
            cfg.OFout = true;
            [cfg,depvars,indepvars]         = SelectSolver(cfg,exvars,depvars,indepvars,r, cfg.CalibrationYears,{'indepvars.ShareWeight'}...
            ,permute(indepvars.ShareWeight(r,1,indepvars.UseModes(r,1,:)),[2 3 1]),cfg.Solver,boundaries.low,boundaries.high);
             end
         end
       
end
end